%clear all; clc;close all;
%computes policy functions for government spending and current account
%Updated: 12/14/2021

label_graph='gn';
label_model='default';

a   =load('agrid2.txt');
b   =load('bgrid2.txt');

vcon0=load('gvcon.txt');
vaut0=load('gvaut.txt');

p0  =load('gp.txt');
h0  =load('gh.txt');
ibp0=load('gibp.txt');
ct0 =load('gct.txt');
ca0 =load('gca.txt');
gn0 =load('ggn.txt');
tau0=load('gtau.txt');
w0  =load('gw.txt');
def0=load('gdef1.txt');
q0  =load('gq.txt');

paut0  =load('gpaut.txt');
haut0  =load('ghaut.txt');
ctaut0 =load('gctaut.txt');
gnaut0 =load('ggnaut.txt');
tauaut0=load('gtauaut.txt');
waut0  =load('gwaut.txt');


na=rows(a);
nk=1;
ns=na/nk;
nb=rows(b);

alpha  = 0.75;
pm     = 1;
gamma  = 0.43;
omegac = 0.3;
mu     = 1;
lambda = 0.184;
r      = 0.02;

atrans = ones(nb,1)*a';

for i=1:na;
    for j=1:nb;
        vcon(i,j)  =vcon0((i-1)*nb + j);
        
        p(i,j)  =p0((i-1)*nb + j);
        h(i,j)  =h0((i-1)*nb + j);
        ibp(i,j)=ibp0((i-1)*nb + j);
        ct(i,j) =ct0((i-1)*nb + j);
        ca(i,j) =ca0((i-1)*nb + j);
        gn(i,j) =gn0((i-1)*nb + j);       
        tau(i,j)=tau0((i-1)*nb + j);       
        w(i,j)  =w0((i-1)*nb + j);   
        def(i,j)=def0((i-1)*nb + j);   
        q(i,j)  =q0((i-1)*nb + j);   
        rec(i,j)=w(i,j)*h(i,j)*tau(i,j);
    end;
end;

bp = b(ibp);

for i=1:na;
    vaut(i)  =vaut0(i);
    
    paut(i)  =paut0(i);
    haut(i)  =haut0(i);
    ctaut(i) =ctaut0(i);
    gnaut(i) =gnaut0(i);       
    tauaut(i)=tauaut0(i);       
    waut(i)  =waut0(i);   
    recaut(i)=waut(i)*haut(i)*tauaut(i);
	wwaut(i) = alpha*paut(i).*haut(i).^(alpha-1); 
    
end;

cn       = h.^alpha-gn;
cnaut    = haut.^alpha-gnaut;
const    = (omegac.*ct.^(-mu)+(1-omegac).*cn.^(-mu)).^(-1/mu);
constaut = (omegac.*ctaut.^(-mu)+(1-omegac).*cnaut.^(-mu)).^(-1/mu);

for i=1:na;
    temp1 = find(ct(i,:) > 0);
    temp2 = find(cn(i,:) > 0);
    temp3 = find(p(i,:) > 0);

    iminct(i) = temp1(1);
    imincn(i) = temp2(1);
    iminp(i) = temp3(1);
    
    imin(i)  = max(iminct(i),max(imincn(i),iminp(i)));
    
    v(i,1:imin(i)-1)  = NaN; 
    p(i,1:imin(i)-1)  = NaN; 
    bp(i,1:imin(i)-1) = NaN; 
    ibp(i,1:imin(i)-1)= NaN; 
    ctE(i,1:imin(i)-1) = NaN; 
    ctU(i,1:imin(i)-1) = NaN; 
    cnE(i,1:imin(i)-1) = NaN; 
    cnU(i,1:imin(i)-1) = NaN; 
    constE(i,1:imin(i)-1) = NaN;     
    constU(i,1:imin(i)-1) = NaN;     
    gn(i,1:imin(i)-1) = NaN; 
    h(i,1:imin(i)-1)  = NaN; 
    ca(i,1:imin(i)-1)  = NaN; 
    tau(i,1:imin(i)-1)= NaN; 
	rec(i,1:imin(i)-1)= NaN;
	ww(i,1:imin(i)-1)= NaN;
	q(i,1:imin(i)-1)  = NaN;
  
end;

for i=1:ns;
for j=1:nk;        
    vcon2(i,j,:)  = vcon((i-1)*nk+j,:); 
    p2(i,j,:)  = p((i-1)*nk+j,:);  
    bp2(i,j,:) = bp((i-1)*nk+j,:);  
    ibp2(i,j,:)= ibp((i-1)*nk+j,:);  
    ct2(i,j,:) = ct((i-1)*nk+j,:);  
    cn2(i,j,:) = cn((i-1)*nk+j,:);  
    const2(i,j,:) = const((i-1)*nk+j,:);         
    gn2(i,j,:) = gn((i-1)*nk+j,:);  
	ca2(i,j,:) = ca((i-1)*nk+j,:);  
    h2(i,j,:)  = h((i-1)*nk+j,:);   
    tau2(i,j,:)= tau((i-1)*nk+j,:);  
	rec2(i,j,:)= rec((i-1)*nk+j,:); 
    def2(i,j,:)= def((i-1)*nk+j,:);  
	q2(i,j,:)  = q((i-1)*nk+j,:); 
	ww2(i,j,:) = alpha*p2(i,j,:).*h2(i,j,:).^(alpha-1); 
    

    vaut2(i,j)  = vaut((i-1)*nk+j); 
    paut2(i,j)  = paut((i-1)*nk+j);  
    ctaut2(i,j) = ctaut((i-1)*nk+j);  
    cnaut2(i,j) = cnaut((i-1)*nk+j);  
    constaut2(i,j) = constaut((i-1)*nk+j);          
    gnaut2(i,j) = gnaut((i-1)*nk+j);  
    haut2(i,j)  = haut((i-1)*nk+j);  
    tauaut2(i,j)= tauaut((i-1)*nk+j);  
	wwaut2(i,j) = alpha*paut2(i,j).*haut2(i,j).^(alpha-1); 
end;
end;

nettax2    = p2.*gn2-ca2;           %net taxes
nettaxaut2 = paut2.*gnaut2;         %net taxes

fiscal2 = tau2 - p2.*gn2;           %fiscal surplus
fiscalaut2=tauaut2-paut2.*gnaut2;   %fiscal surplus

nfiscal2 = nettax2 - p2.*gn2;           %fiscal surplus net of transfers
nfiscalaut2=nettaxaut2-paut2.*gnaut2;   %fiscal surplus net of transfers

lowa = round(2*ns/6);
meda = round(ns/2);
hgha = round(4*ns/6);

b    =-b/(lambda+r);  %change of sign
bdef = b;
a    = exp(a);

%put together policy fcns of repayment and autarky
mat_ones = ones(size(def));
p3  =def.*p+(mat_ones-def).*(paut'*ones(1,nb));
h3  =def.*h+(mat_ones-def).*(haut'*ones(1,nb));
ct3=def.*ct+(mat_ones-def).*(ctaut'*ones(1,nb));
cn3=def.*cn+(mat_ones-def).*(cnaut'*ones(1,nb));
ca3 =def.*ca+(mat_ones-def).*0;
gn3 =def.*gn+(mat_ones-def).*(gnaut'*ones(1,nb));       
tau3=def.*tau+(mat_ones-def).*(tauaut'*ones(1,nb));       
w3  =def.*w+(mat_ones-def).*(waut'*ones(1,nb));   
rec3=def.*rec+(mat_ones-def).*(recaut'*ones(1,nb));

[i1,i2] = find(def>0);
q3  =NaN; bp3 =NaN; sp3 =NaN; gn3 = NaN; h3=NaN; p3=NaN;

for i=1:na;
for j=1:nb;
    if (def(i,j)>0); q3(i,j)=q(i,j);bp3(i,j)=bp(i,j);sp3(i,j)=100*((1+q3(i,j).^(-1)-lambda)./(1+r)-1);p3=p;h3=h;gn3=gn;tau3=tau;w3=w;rec3=rec;
    else
        q3(i,j)=NaN;bp3(i,j)=NaN;sp3(i,j)=NaN;p3=NaN;h3=NaN;gn3(i,j)=gnaut(i);tau3=NaN;w3=NaN;rec3=NaN;
    end;
end;
end;

nettax3    = p3.*gn3-tau3;          %net taxes
fiscal3    = tau3 - p3.*gn3;        %fiscal surplus
nfiscal3   = nettax3 - p3.*gn3;     %fiscal surplus net of transfers


NameArray = {'LineStyle','Color','LineWidth'};
NameArray0 = {'LineStyle','Color','LineWidth'};
ValueArray0 = {'--','-';[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5}';
ValueArray1 = {'-';[1  0.4  0.4];1.5}';
ValueArray = {'--','-';[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5}';
ValueArray2 = {'--','-','--','-';[1  0.4  0.4],[0.078  0.1686  0.549],[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5,1.5,1.5}';
ValueArray4 = {'--','--','-','-';'red','blue','red','blue';1.5,1.5,1.5,1.5}';
ValueArray5 = {':',':';[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5}';

for i=1:na;
    for j=1:nb;
        if (def(i,j)==0); 
            ca3(i,j)  =NaN;
            gn3(i,j)  =gnaut(i);
        end;
    end;
end;

load_rf_fcns
cd ..
cd ..
cd ..

endb = 2;
lowb = 159;
hghb = 122;    %average debt: index=122 value=-1.108
hghb2 = 1888;

NameArray = {'LineStyle','Color','LineWidth'};
ValueArray4 = {'--','-','--','-',;[0.078  0.1686  0.549],[0.078  0.1686  0.549],[0.8  0  0],[0.8  0  0];2,2,2,2}';  
ValueArray1 = {'--';[0.078  0.1686  0.549];2}';
ValueArray = {'--','-',;[0.078  0.1686  0.549],[0.8  0  0];2,2}';  

%plot government spending as function of debt
gn3(9,77)=NaN;


%plot government spending as function of yT
gn4(1:4)=gn3(1:4,hghb);
gn4(5)=0.21935;
gn4(6)=NaN;
gn4(7:na+2)=gn3(5:na,hghb);
a4(1:4)=a(1:4);
a4(5)=a(5)-0.0001;
a4(6)=a(5)-0.001*0.01;
a4(7:na+2)=a(5:na);

nameFig = char([char(label_graph) '_yt']);
figure('name',nameFig);
P = plot(a4,squeeze(gn4),a,squeeze(gnrf(:,hghb2)));
lgd=legend('defaultable debt','risk-free debt',' ','FontSize',19);lgd.Location='SouthEast';legend boxoff;
hold on;plot([a(5) a(5)],[0 1],':k','LineWidth',2);
ylabel('g^{N}','Fontweight','normal','FontSize',16); 
xlabel('y^T','Fontweight','normal','FontSize',16);
AX1 = gca;box on;
set(AX1,'YTick',0.15:0.05:0.25,'FontSize',16);
axis(AX1, [a(1) a(end) 0.15 0.25]);
set(P,NameArray,ValueArray);grid on;
saveas(gcf, fullfile(plot_path, nameFig),format_chart);
